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Abstract 

In this work, a complete first order analytical solution, which includes the short periodic terms, for the 
problem of optimal low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian 
central gravity field is obtained through Hamilton- Jacobi theory and a perturbation method based on Lie series. 

1. INTRODUCTION 

In the general analysis of optimal space trajectories, two idealized propulsion models are of most 
frequent use [1]: LP and CEV systems. In the power-limited variable ejection velocity systems - LP systems - 
the only constraint concerns the power, that is, there exists an upper constant limit for the power. In the constant 
ejection velocity limited thrust systems - CEV systems - the magnitude of the thrust acceleration is bounded. In 
both cases, it is usually assumed that the thrust direction is unconstrained. The utility of these idealized models is 
that the results obtained from them provide good insight about more realistic problems. 

The purpose of this work is to present a complete first order analytical solution, which includes the 
short periodic terms, for the problem of optimal low-thrust limited-power simple transfers (no rendez-vous) 
between arbitrary elliptic coplanar orbits in a Newtonian central gravity field. This analysis has been motivated 
by the renewed interest in the use of low-thrust propulsion systems in space missions verified in the last two 
decades. 

The optimization problem is formulated as a Mayer problem of optimal control theory with Cartesian 
elements - position and velocity vectors - as state variables [1]. After applying the Pontryagin Maximum 
Principle and determining the maximum Hamiltonian, classical orbital elements are introduced through a 
canonical transformation - Mathieu transformation. The short periodic terms are then eliminated from the 
maximum Hamiltonian through an infinitesimal canonical transformation built through Hori method - a 
perturbation canonical method based on Lie series [2]. The new maximum Hamiltonian, resulting from the 
infinitesimal canonical transformation, describes the extremal trajectories associated with the long duration 
maneuvers. Using Hamilton- Jacobi theory [3], closed-form analytical solutions are obtained for the canonical 
system governed by this new Hamiltonian. For long duration maneuvers, the existence of conjugate points is 
investigated through the Jacobi condition and the envelope of extremals is obtained considering different 
configurations of initial and final orbits. Curves representing maneuvers with same fuel consumption for a 
specified duration, called iso-cost curves, are presented for several configurations of initial and final orbits. The 
results are very similar to the ones obtained by Edelbaum [4] and Marec and Vinh [5]; but, the theory presented 
in the paper includes directly the periodic terms such that it can be applied to any transfers independently of its 
duration. On the other hand, since the theory is based on canonical transformations, a first order solution for the 
Cartesian elements can also be easily obtained from Lie theorem [2]. 

2. OPTIMAL SPACE TRAJECTORIES 

Let us to consider the motion of a space vehicle M, powered by a limited-power engine in a Newtonian 
central gravity field. At time t, the state of the vehicle is defined by the position vector r(t) , the velocity vector 
v{t) and the consumption variable J [1], 

j^U'fdi, (i) 

2 J l <> 


i 



where y denotes the magnitude of the thrust acceleration y , used as control variable. It is assumed that the 
thrust direction and magnitude are unconstrained. 

The optimization problem will be formulated as a Mayer problem of optimal control as follows. It is 
proposed to transfer the space vehicle M from the initial state (r 0 , v 0 ,0) at the initial time t 0 = 0 to the final state 
( r, , v , , J , ) at the specified final time i , such that the final consumption variable J . is a minimum. The state 
equations are 
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where p is the Gaussian constant. 

According to the Pontryagin Maximum Principle [6], the optimal thrust acceleration y must be 
selected from the admissible controls such that the Hamiltonian H reaches its maximum. The Hamiltonian is 
formed using Eqns (2), 


H = p r • v + p 
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where p r , p and p , are the adjoint variables and dot denotes the scalar vector product. Since the optimization 
problem is unconstrained, it follows that 

y'=~PjPj • ( 4 ) 


The optimal thrust acceleration j i is modulated [1], Then, the optimal trajectories are governed by the 
maximum Hamiltonian H‘ , 
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Since J is an ignorable variable, it is determined by simple integration, and its adjoint Pj is a first integral, 
whose value is obtained from the transversality conditions: p , = - 1 . Thus, 
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H' = p r .v- p v .-^r + ^~. (5) 
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This Hamiltonian can be decomposed into two parts: H 0 = p r • v - p v • r is the undisturbed Hamiltonian and 
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H . = — is the disturbing function concerning the optimal thrust acceleration. 

3. TRANSFORMATION FROM CARTESIAN ELEMENTS TO ORBITAL ELEMENTS 

In order to obtain a first order analytical solution for simple transfers between arbitrary elliptic coplanar 
orbits, the set of classical orbital elements is introduced through a canonical transformation by using the 
properties of generalized canonical systems [7]. Consider the extended two-body problem defined by the 
canonical system of differential equations governed by the undisturbed Hamiltonian H 0 , that is, the problem 
defined by the differential equations for the state variables r and v , which describe the classical two-body 
problem, and the differential equations for the adjoint variables p r and p , 

dt dt r 
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where e r is the unit vector pointing radially outward of the moving frame of reference. In the coplanar case, the 
general solution of this system is given by [7]: 
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where e i is the unit vector along circumferential direction of the moving frame of reference; a is the semi-major 
axis, e is the eccentricity, a> is the pericenter argument, / is the true anomaly, E is the eccentric anomaly, M is 
the mean anomaly, n = -Jp/ci 3 is the mean motion, and ( r/a ), (r/a)sin/, ... etc are known functions of the 

elliptic motion, which can be written explicitly in terms of eccentricity and mean anomaly through Lagrange 
series [8]. The anomalies are related through the equations: 

f ll + e E 

tan— = tan — , and, M = E -esin£ (Kepler's equation). 
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The unit vectors e r and e s are written as function of the orbital elements as follows 

e r = [cosf2cos(o + /)-sinQsin(ft) + /)cos/]i +[sinQcos(® + /) + cosf2sin(o + /)cos/}/ 

e s = -[cosf2sin(® + /) + sinQcos(® + /)cos/]i + [-sinQsin(® + /) + cosflcos(® + /)cos/]y 

where i and j are the unit vectors of a fixed frame of reference. 

Equations (6) - (9) define a canonical transformation - Mathieu transformation - between the Cartesian 
elements (r,v, p r , p r ) and the orbital ones (a,e,co,M ,p a ,p e ,p (o ,p M ) . The Hamiltonian El * is invariant with 
respect to this canonical transformation. Thus, 

H 0 =np M , (10) 
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The new Hamiltonian function //*, defined by Eqns (10) - (11), describes the optimal low-thrust 
limited-power trajectories in a Newtonian central gravity field for transfers between arbitrary elliptic coplanar 
orbits. 


4. A FIRST ORDER ANALYTICAL SOLUTION 


In order to get an approximate formal solution for the problem of transfers between arbitrary elliptic 
coplanar orbits, a perturbation technique based on Lie series - Hori method - will be applied. For completeness, 
a brief description of Hori method [2] is presented in the Appendix. 

In this study, it is assumed that the functions // and H . have different order of magnitude for low- 

thrust propulsion systems. So, from Eqns (10) - (11), using the algorithm described in the Appendix, one finds 
the generating function S and the new Hamiltonian F ' , given up to the first order by 
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Prime denotes the new canonical variables resulting from the infinitesimal canonical transformation generated by 
S. Terms factored by p' M do not appear in equations above, since only simple transfers are considered [4]. It 
should be noted that mixed secular terms can arise in the new Hamiltonian in the case of rendez-vous maneuvers. 
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The new Hamiltonian F' governs the optimal long duration transfers between arbitrary elliptic 
coplanar orbits. The general solution of the canonical system described by this Hamiltonian will be obtained 
through two canonical transformations. First, consider the Mathieu transformation defined by the following 
equations: 
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The Hamiltonian function F' is invariant with respect to this transformation. Thus, 
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The canonical system governed by the new Hamiltonian function F" has three first integrals (t and m 
are ignorable variables), 
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These first integrals play important role in determining the general solution of the canonical system governed by 
the new Hamiltonian function F" through Hamilton- Jacobi theory. 

Let us to consider the problem of determining a canonical transformation defined by a generating 
function W such that C, , C 2 and E become the new generalized coordinates, 


(a", <j>, a”; p”, p* , p " a ) > (c, ,C 2 ,E;p c ,p^,p E ). 


Since the new Hamiltonian function F" is a quadratic form in the conjugate momenta, the separation of 
variables technique will be applied for solving the Hamilton- Jacobi equation [3]. 

Consider the transformation equations 
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where W = W(a",^,co",C 2 ,C 2 ,E) is the generating function. The corresponding Hamilton- Jacobi equation is 
then given by 
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Following the separation of variables method, it is assumed that the generating function W is equal to a sum of 
functions, each of which depends on a single variable, that is: 

W(a",0,a)",C 1 ,C 2 ,E) = W l (a’,C l ,C 2 ,E) + W 2 (<f,C l ,C 2 ,E) + W 3 (a>’,C l ,C 2 ,E) . (21) 

Therefore, from Eqns (16) through (21), it follows that 
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A complete solution of these equations is given by 
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W 2 =j yjc 2 2 - C; esc 2 <t d </> , (23) 

W 3 = C.ffl" , (24) 


with 5C 2 -4Cj 2 = 5C 2 . We note that W 2 is given as indefinite integral, since only its partial derivatives are 
needed (Eqns (19)), as shown in the next paragraphs. 

Now, consider the differential equations for the conjugate momenta of the canonical system governed 
by the new Hamiltonian F'" = E , that is: 

d Pc t 0 d Pc 2 Q dp E 

dt dt dt 

whose solution is very simple: 

Pc 1 — ^ Pc 2 — Pe “ ^3 — ^ 5 (25) 

where a„i = 1,2,3 , are arbitrary constants of integration. 

Introducing the generating function W, defined by Eqns (21) through (24), into the transformation 
equations, Eqns (19), and taking into account the general solution defined by Eqns (25) for the conjugate 
momenta, one gets 
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From equations above, one finds the solution of the canonical system governed by the Hamiltonian F" 
for a given set of initial conditions: 
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with the auxiliary constants k 0 , k t and k 2 defined as functions of the initial value of the adjoint variables 
(conjugate momenta) by 


esc 2 k„ = 


s^oX) 2 +X"( 5csc2/t , - 4 ) 

X( 5csc;2 x 4 ) 


esc 2 k, = 


Pl+P^csc 2 ^ 


k 2 =®''+— r 0 sin k t - tan 1 (tan r 0 esc k t ) . 


The constants C, Cj , C, and E in Eqns (26) through (31) can also be written as functions of the initial value of 
the adjoint variables as follows: 

C 2 =^K;(5csc 2 A: 1 -4), C, = p\ , 

c i= Pl+P'^ csc2 A’ 4pE = a" 0 (p (a" 0 p" ao ) + p” 2 (5csc 2 k t -4)). 

The initial conditions for the state variables (generalized coordinates) are given by a"(o)= a " , e(0) = sin and 
®"(o) = a >" , and, r 0 is obtained from cos(*i 0 = cos A, cosr 0 . 

Equations (32) through (39) represent the solution of the canonical system concerning the problem of 
optimal long duration low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a 
Newtonian central field. 

For maneuvers with arbitrary duration, the periodic terms must be included [9]. Following Hori method, 
a first order analytical solution is then given by: 
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with a’ , e , p’ m given through Eqns (14). Equations (40) - (42) are in agreement to the ones obtained in [9] 
through Bogoliubov-Mitropolsky method (a non-canonical perturbation technique). The eccentric anomaly E' is 
computed from Kepler's equation with the mean anomaly M' given by 
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where M' 0 =M'(t 0 ). The differential equation for the mean anomaly M' is derived from the undisturbed 
Hamiltonian F 0 and an additional term A F' factored by p' M , 
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where dots denote the terms defined by Eqn (12). 

For transfers between coaxial orbits with no change in the argument of pericenter, the analytical 
solution described in the preceding paragraphs simplifies. Since p’ = 0 , some equations must be rewritten: 
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5. ANALYSIS OF LONG DURATION TRANSFERS 

In this section, a complete analysis of long duration transfers described by the Hamiltonian F' is 
presented. The investigation of conjugate points and the solution of the two-point boundary value problem are 
considered. 
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To perform this analysis, a new consumption variable is introduced [4], u = ^2E t , and Eqns (32) and 
(33) can be put in the form: 
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where v 0 = yJ/j/a’ . By eliminating k 0 from these equations, one finds: 



On the other hand, J = El . thus: 
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Now, consider Eqns (33) - (36) which define a two-parameter family of extremals in the phase space 
(fl", <j>, a") for a given initial phase point (a",^ 0 ,ai") corresponding to an initial orbit. By eliminating the auxiliary 
variables t andy/ , a = a" /a" and a" can be written as explicit functions of (f> , that is, a = a(<j>,</> 0 ,k 0 ,k^) and 
a" = co"{<l),<j> 0 ,(o i) ,k^) . The conjugate points to the phase point (a", $,,©') are determined through the roots of the 
equation [10]: 
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analysis of the roots of the following equation: 


from Eqns (33) and (47), one finds 
the conjugate points reduces to the 



(50) 


From Eqns (35) and (36), one finds the explicit form of a" = a"(if>,<j> 0 ,a 0 ,k l ) , that is, 
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The partial derivative da"/dk x is given by 
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The auxiliary variable r is re-introduced to simplify the expression. Therefore, from Eqns (50) and (51), after 
some simplifications, it follows that the conjugate points are given by the roots of following equation [4], 
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As discussed in [11,12], the conjugate points determined through Eqn (52) occur for transfers during which the 
direction of motion in the orbit reverses; that is, for transfers between direct orbits (no direction reversals) the 
extremals computed for the canonical system governed by the Hamiltonian F' are globally minimizing. 

To complete the analysis, a brief description of an iterative algorithm for solving the two-point 
boundary value problem of going from an initial orbit O 0 : ( a",e",a ?”) to a final orbit O , : (a " , e " , co” ) is presented. 
The steps of this algorithm can be described as: 


1 . 

2. 

3. 

4. 

5. 

6. 
7. 


Guess a starting value of k , ; 

Determine (f n , </ f , r 0 and r, through Eqns (14) and (35); 
Compute co" through Eqn (36). 


If co" * co " f , adjust the value of k ] and repeat steps 2 and 3 until \co” -©*<£, where s is a prescribed 
small quantity; 

Compute y/ f and u f /v 0 through Eqns (34) and (48); 
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The adjust of the value of k x can be obtained through the well-known Newton-Raphson method, that is, 
co"{k l „)-£», 

{dcD'/dkX^ ’ 


with the partial derivative (dco"/dk l \^ computed from Eqn (51). In the case of coaxial orbits, k, = 0 and the 

solution of the two-point boundary value problem involves only the solution of Eqn (43), which is given through 

step 6 with replaced by — j=((zS-(zS 0 ) . The initial value of the adjoint variables is given as follows: p’ is 

•v 5 

computed through step 7 and p h is given by p 1 ^ = 

The thrust acceleration required to perform the maneuver is given through Eqns (4) and (9), that is, 
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with the state and adjoint variables computed from Eqns (32) through (39). 
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6. RESULTS 


In this section some numerical results are presented through Figs 1-5. Figures 1 and 2 show the field 
of extremals and iso-cost curves for long duration transfers between coaxial orbits, considering two cases: in the 
first one, the space vehicle departs from a circular orbit and, in the second one, it departs from an elliptic orbit 
with eccentricity e 0 = 0.4 . In both cases, the semi-major axis is a 0 = 1 , such that the results are presented in 
canonical units. The extremals are plotted for different values of the constant k 0 and the iso-cost curves for 
different values of the new consumption variable u/v 0 . 

120 ° 90 = 60 " 45 " 30 " 15 " 5 " 



Figure 1 - Field of extremals and iso-cost curves for transfers between coaxial orbits for e 0 = 0.0 . 



Figure 2 - Field of extremals and iso-cost curves for transfers between coaxial orbits for e 0 = 0.4 . 
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Figures 3 and 4 show the field of extremals and iso-cost curves for long duration transfers between non- 
coaxial orbits. Note that in the general case the family of extremals is described by two constants k 0 and /c . In 
the plots of Figure 3, k 0 = 18° has been chosen in order to illustrate the general aspect of the field of extremals 
and compare the results with the ones obtained in [4,5,11,12]- The projection of extremals on the (e,®) -plane 
shows clearly the non-existence of conjugate points for the extremals computed from Eqns (32) - (36) in the 
case of transfers between direct orbits (no direction reversals). Figure 4 shows a sample of iso-cost curves for 
different values of the new consumption variable w/v 0 and a fixed value of = 18° . Note that these curves lie 
on an iso-cost surface with complex geometry. By simplicity, the symmetric interval [-180°, 180°] for co is 
considered. 




Figure 3 - Field of extremals for transfers between non-coaxial orbits for e 0 = V 3/2 with k 0 =18° 



Figure 4 - Iso-cost curves for a 0 =1.0, e 0 =0.5 and ® 0 = 0° , with k x =18°. 

Figure 5 shows a comparison between the complete analytical solution, the average solution and the one 
obtained through numerical integration of the set of canonical equations which describes the optimal trajectories 
(system of differential equations governed by new Hamiltonian function H\ defined by Eqns (10) - (11)) for a 
transfer between coaxial orbits with the departure orbit O 0 defined by the set of initial conditions a 0 =1.0, 
e 0 = 0.20 , ® 0 = 0° and the arrival orbit defined by the set of terminal conditions a 0 = 2.0 , e 0 = 0.25 , a> 0 = 0° , 
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and a transfer duration t f —t 0 = 500 canonical units. The two-point boundary value problem is solved through a 

Newton-Raphson algorithm based on the complete analytical solution presented in Section 4 which uses the 
algorithm described in the preceding section for long duration transfers to generate a starting approximation of 
the adjoint variables. Note that there exists a good agreement between the complete analytical solution and the 
numerical one; but, the average solution loses accuracy for large values of the time when compared to the others 
solutions. 

Finally, we note that a first order solution for the Cartesian elements can also be easily obtained from 
Lie theorem [2], since the theory presented in the paper is based on canonical transformations. 


6. CONCLUSION 

A first order analytical solution for the problem of optimal low-thrust limited-power transfers between 
arbitrary elliptic coplanar orbits in a Newtonian central gravity field has been obtained through Hamilton-Jacobi 
theory and a perturbation method based on Lie series. A complete analysis of long duration transfers, which 
includes the investigation of conjugate points and the solution of the two-point boundary value problem, is 
presented. Results show the good agreement between the complete analytical solution and the one obtained 
through numerical integration of the set of canonical equations which describes the optimal trajectories for a 
transfer between coaxial orbits. On the other hand, the average solution loses accuracy for large values of the 
time when compared to the complete analytical solution. 





1 


2.4 


Figure 5 - Comparison between numerical, analytical and average solutions for a transfer between coaxial orbits. 
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APPENDIX 

Hori method is based on a Lie theorem [2]. 

Let H{x,y,s ) be a Hamiltonian written as a power series of a small parameter s. 
H(x,y,s) = H 0 (x,y) + Y J £ t H t (x,y) , 


where the undisturbed Hamiltonian H 0 (x,y) describes an integrable system. Here, x and y are n x 1 -vectors of 
generalized coordinates and conjugate momenta, respectively. 

Consider an infinitesimal canonical transformation (x, y) — > (%,rj) defined by a generating 

function S(g, t], s) , also developed as a power series of the small parameter s , 

k = 1 

This transformation is such that the new dynamical system has some advantages for the solution, f and ij are 
also n x 1-vectors. The new Hamiltonian F(d.ij.c) , resulting from this canonical transformation, is also 
developed as a power series of the small parameter s , 

F^,i],£)=F 0 ^,i]) + Y J £ k F t ^,'l)- 

k=\ 


The transformation equations are given explicitly by 



where D,."f = f, D s 'f = {/,S}, D"f = D" 1 (a. 1 /) > n>2, f{§,T]) denotes an arbitrary function and braces 
stand for Poisson brackets. 

According to the algorithm of Hori method, the new Hamiltonian F(^,rj,s) and the generating 
function S(c.)].c) are obtained, at each order of the small parameters , from the equations: 
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• order 0: H a (£ ,ij) = F a (£ ,rj) , (A.l) 

• order 1: {H 0 , S, }+//,= F, , (A.2) 

• order 2: +F X ,S)+H 2 = F 2 . (A.3) 

Functions in above equations are written in terms of the new set of canonical variables (c,//) . 

The »;-th order equation of the algorithm can be put in the form 

= (a.4) 


where 0 m is a function obtained from the preceding orders, involving H 0 , H m , .S', , F k and H k , k = \,...,m-\. 

The determination of the functions F m and S m is based on the general solution of the undisturbed 

system 


dt dij. dt 


(A.5) 


and it is performed through an averaging principle. Following the integration theory proposed in [13], Eqn (A.4) 
can be put in the form 


dS m 

dt 


= ©,„ 


- F 


with 0 m written in terms of the arbitrary constants of integration, c 1 ,...,c u , of the general solution of the 
canonical system (A.5). S m and F m are unknown functions. 

Following Flori [2], suppose that the canonical transformation generated by the function S is such that 
the time t is eliminated from the new Flamiltonian F. This is accomplished taking F m as the mean value of the 
function @ m . Therefore, functions S m and F m are given through the equations 


F„=(® m ), (A. 6) 

S m =J [©.-(©„ ))lt, (A. 7) 

where ( ) stands for the mean value of the function. It should be noted that the averaging process and the 
integration in Eqns (A.6) and (A.7) are performed considering the explicit time dependence of the functions. 
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